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1. Introduction 

RNA molecules play an important role in all living organisms [Ij. They are usually 
found in a at least partially folded state, due to the pairing of a base with at most one 
other base. A given configuration is thus characterised by the set of base pairings, see 
figure [H These pairings are mostly planar [21 |3l H] (see [5] for non-planar corrections), 
which is what we will suppose from now on. At high temperatures, in the so-called 
"molten phase", energetic considerations only play a minor role, and the probability Pij 
of two RNA-bases to pair, is [6] 

where i and j are the labels of the bases counted along the backbone/strand, and n is 
the overall size of the RNA-molecule, i.e. its total number of bases. 

At low temperature, the RNA-molecule will settle into the optimally paired (or 
folded) configuration, i.e. the minimal energy state, as long as this state is reachable 
in the available time-scales. The optimal fold for a given molecule is a question to be 
answered by biology. Since biological sequences are rather specific, much effort has been 
invested to understand the properties of a random sequence, termed "random RNA". 
The idea is that either the folding properties of random RNA are close to those of 
biological sequences, or if not, that they must be characterised in order to understand the 
deviations present for biological RNA, giving eventually a hint why nature is organising 
in a certain way. 




(a) (b) (c) 

Figure 1. RNA molecules, like DNA, are long chain heteropolymers built from 
four types of nucleotides: adenine (A), uracil (U), guanine (G) and cytosine (C). In 
solution, a single RNA molecule bents back onto itself and folds into a configuration 
of loops, stems and terminating bonds, due to pair formation from nucleotides located 
on different parts of the polymer strand. The set of base pairs, Watson-Crick pairs 
A-U, G-C and the less favorable wobble pair G-U, defines the secondary structure. 
Illustration of RNA secondary structures: (a) an RNA molecule with given base 
sequence folds into a base pair configuration (b). In the absence of pseudo-knots 
the secondary structure may be represented as a diagram of non-intersecting arches 
(c). 
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Characterising random RNA has proven a challenge so far: Numerical work 
[21 m [HI [9] is restricted to relatively small molecules, with up to maximally 2000 bases 
[10] . despite the fact that rather efficient polynomial algorithms exist (~ n^). Analytical 
work was pioneered by Bundschuh and Hwa [21 [7]. From their numerical work, they 
claim that for large molecules, a random-base model is indistinguishable from a random 
pairing-energy model, where the pairing energy e^j between base i and j is a random 
Gaussian variable, confirmed in [8]. Bundschuh and Hwa then conjectured that a phase 
transition separates the high-temperature molten phase from a low-temperature frozen 
phase. Using an RG treatment, Lassig and Wiese [11] showed analytically, that this 
phase transition exists, and is of second order. They also calculated the exponents 
characterising the transition, and using a locking argument extended their findings to 
the low-temperature (glass) phase. David and Wiese [T2] substantiated these findings, by 
constructing the field theory, showing its renormalisability to all orders, and performing 
an explicit 2-loop calculation, yielding 

Ptransition Pfrozen ~ 1.36 . (2) 

The field theory makes some definite predictions about the transition, which are hard to 
verify numerically. A major problem is that the systems are not large enough to analyze 
the asymptotic behavior. Under these circumstances, the knowledge of a scaling function 
would be very helpful, as would be the knowledge of the form of corrections to scaling. 

We therefore propose a simple hierarchical model, where all this can be calculated 
analyticalljljl. This is based on the observation, that if the n{n — l)/2 possible pairing 
energies ejj are ordered hierarchically 

then the construction of the minimal energy configuration is much simplified: First take 
the largest pairing energy e^jj^, and pair bases ii and ji. Among the remaining pairings, 
consider only those allowed by planarity. Among those, choose the one with the largest 
pairing energy, and pair the corresponding bases. Repeat this procedure until no more 
bases can be found. The same idea is at the base of the dynamics for greedy algorithms 
of RNA folding: At each time-step, choose the most favorable base pairing and fold it. 

In this article, we systematically analyse the statistical properties of the structures 
built in the hierarchical model. In particular, we compute exactly its properties as 
n — > oo, for example we prove that in this limit the pairing probability reads 

^(^,j)~N-jr^ P = ^ - 1-44 . (4) 

We then calculate scaling functions for higher moments of the "height function" 
(which encodes the pairings), and their finite-size corrections. This is achieved with 
two complimentary approaches: Generating functions for the arch-deposition model 

I While working on this project, we learned from Markus Miiller that he had considered this model in 
his PhD-thesis [13], but not published elsewhere. He also foimd the scaling exponent ( to be discussed 
below, but did not consider the scaling-functions and corrections-to-scaling which are the main purpose 
of this article. 
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introduced above, and a dual tree-growth process. The advantage is that quantities 
which can easily be calculated in one model, are difficult to obtain in the other, and vice 
versa. This idea may be interesting for more general tree-growth processes, since if the 
dual model can be constructed there, it would allow to calculate otherwise inaccessible 
quantities. For examples of tree growth processes, we refer the reader to [m [15], HH [17] 
among the vast existing literature. 

The presentation is organised as follows. In section [2], we provide a general 
framework for the hierarchical model in terms of recursion relations for finite n. The 
recursion relations are analysed by means of generating functions. In the limit n oo 
we extract the scaling behaviour of various quantities and compute sub-leading finite- 
size corrections in sections [3], H] and [5l We compare our results to numerical simulations 
in section [6] In section [7] we present an alternative tree-growth model which we show 
to yield equivalent structures even though the dynamics of their construction is quite 
different. Several technical points and extensions are relegated to three appendices. 



2. Arch deposition model 

2.1. Arch systems and height functions 

We consider a strand with n bases labeled by indices i = 1, . . . , n. Similarly, we use the 
same index i to label the segments between consecutive bases i and i + 1. A secondary 
structure C is a set of base pairs with 1 < i < j < n. C is called planar if any 
two {k,l) G C are either independent i<j<k<loT nested i < k < I < j. In 

what follows, the structures are supposed to be planar. Thus, we may represent a given 
structure by a diagram of non-intersecting arches (see figure [2^). 




Figure 2. (a) Arch diagram for a planar structure, (b) Corresponding tieight relief, 
defined in 



Given some structure C, it is natural to ask whether it contains an arch a = {i,j). This 
is answered by the contact operator defined by 

^ U j\ — J " if a e C 
' otherwise 
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For our investigations the so-called height function for the segment i will play a central 
role. It is defined as 

i n 

hc{i.n):=Y, E "^cij.k). (6) 

j = l k=i+l 

It counts the number of arches above a given segment [z, i + 1], and thus has boundary 
conditions hc{0,n) = hc{n,n) = 0. Therefore, the height function hc{i,n) provides 
a one-to-one correspondence between C and mountain reliefs (Dyck-like paths) on the 
interval [0,n] subject to vanishing boundary conditions and \hc{i + l,n) — hc{i,n)\ = 1 
or (see figure [2)d). We define the average height by 

^ n—l 

hc{n) := - ^hc{i,n) . (7) 

i=l 

2.2. The random arch deposition process 

Definition of the model (model A). The structures C are built up in the following way: 
At initial time step t = 0, we start with n unoccupied points on the line. At each time 
step t, we deposit a new arch as follows. At time step t — 1, we have already a planar 
system of t — 1 arches linking 2(t — 1) points. We have m = n — 2{t — 1) free points left, 
and we may build m{m — l)/2 different arches. We now consider the subsets of these 
arches {i,j) which keep the arch system planar, when added to the present structure. 
We choose at random, and with equal probability, one of these arches, and add it to the 
system at time t (as depicted on figure [3]). 

The process is stopped as soon as no more planar deposition is possible. The 
stopping time (tstop = number of arches of the final configuration) will vary from 
configuration to configuration, since not all points get paired. 

We call this arch deposition process "model A" . 
• •••••• /^ • ' • ^ ^• \ • m 

Figure 3. Building up planar structures via successive arch deposition. 

Hierarchy and recursion for probabilities. Our construction is "hierarchical" in the 
sense that each deposition partitions the strand into two non-connected substrands. 
Since this procedure is performed at random, it naturally induces a probability measure 
Pa{C) on the set of structures C with a given number of points n. Although it turns out 
to be a quite tedious exercise to compute the probabilities for structures, even with only 
= 4, 5, 6, . . . points, we can write a formal yet powerful recursion relation for Pa{C). 
Given C, any arch a e C may have been the first arch in the construction process (at 
t = 1). Since the deposition of a = («, j), 1 < i < j < n, is not constrained by the 
presence of any other arches, its probability is uniform and simply given by 2/[n(n — 1)]. 
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This first step leads to a separation of tlie strand into an "interior" part witli ni = j—i—2 
points and an "exterior" part witli n2 = n — j + i. Tlie deposition process tlien grows 
structures Ci inside and C2 outside tlie first arch (see figure Hj). The key observation is 
that these structures grow independently. Therefore, their joint probability factorises: 




Figure 4. Decomposition of a configuration C in model A 

This recursion relation, together with the initial condition P4 = 1 for the n = and 
the n = 1 configurations (no point and a single free point), is sufficient to obtain all 
probabilities. In fact, it is this relation that renders the arch-deposition model amenable 
to exact analytic calculations. 

With the help of Pa{C) we can compute averages, i.e. expectation values of an 
observable J-'c via 

•^=(-^c) = $^Pa(C)^c, (9) 

c 

where the sum is carried out over all possible structures with a fixed number n of 
points. Throughout this article, we follow the convention to note objects depending on 
an individual structure C with a subscript. 



2. 3. Summary of results 

In this article, we focus on the mean height at a given point h{i,n) = {hc{i,n)) and 
the probability P{i,j) that two bases located at i and j are paired. The latter is the 
expectation value of the contact operator ^c{hj), P{hj) '■= {^c{hj))- 

Before embarking into calculations, let us briefly summarise some important 
properties of these quantities as well as our main results. First, note that the 
construction of the height function ([6]) implies that h{l,n) is the probability that point 
2 = 1 is paired to any other point 2 < j < n on the strand. Since averaging over 
structures will lead to translational invariance of {^c{hj)) = {^c{i + + ^)), for all 
m, we can interpret h{l,n) as the probability that some arbitrary point 1 < z < n is 
involved in a pair. We compute h{l,n) for any n and show that it converges to 

lim /i(l,n) = 1 -e"^ = 0.864665... (10) 
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The full information about all possible structures is contained both in the height profiles 
as in the pairing probabilities. In the scaling limit n oo, we show that the height 
function and the pairing probabilities take the scaling forms 

h{i,n) ~ n'^Hii-] and P{i,j) ~ n"" P ( ) (11) 

n—>oo V'^/ n-^oo V / 

with scaling functions TCi and V which we compute exactly as well as the scaling 
exponents ( and p. From eqn. (Q, we immediatly deduce the scaling relation ( + p = 2. 

The exponent ( is also related to the intrinsic Hausdorff dimension of the tree 
structure dual to the arch system by dh = 1/C- Therefore it is sufficient to determine 
the exponent ( which we show to be 

C = ^^-^ ^ 0.561553 , p = ^ — ^ 1.43845 . (12) 

This agrees with [13]. It follows that the average mean height h{n) = {hc{n)) grows like 
for large n. We determine its exact generating function, which allows us to compute 
sub-leading corrections to the scaling limit to any desired accuracy. 

The analysis of higher moments {h{i,n)^) naturally raises the question of 
multifractality of the arch structures/height profiles. We show that 

{hc{i,nf) ~ n^mJ-\ withCfc = A;C (13) 

with scaling functions Tik that we can in principle compute. Using this result we are 
able to prove the absence of multifractality. 

Recurrence relations and generating functions 

We now exploit the recursion relation ([8]) to compute the moments of the height function 
{hcii, nY), /c = 0, 1, 2, . . . Our general strategy is to extract their properties by analyzing 
the behavior of their corresponding generating functions. 

2.4-1 ■ Recurrence relation for the height function: the principle We want to evaluate 
the height hc{i., n) for a given structure C. The first arch a = (j, k) splits C into the two 
independent substructures Ci and C2 with lengths ni = n — k + j — 1 and n2 = k — j — 1 
respectively. We now consider the height over segment [i,i + 1]. With respect to first 
arch a, this segment may have three different locations, as indicated on figure [5l (a) if 
i < j, the segment is situated on the part of the strand which belongs to Ci and thus the 
height is given by hc^{i,n — k + j — 1). (b) The case i > k is similar, but we must shift 
the position i — > z — /c + j — 1, we thus find the height hc^{i — k + j — l,n — k + j — 1). 
(c) Finally, if j < i < k, we have to count the height for the structure C2 with the 
readjusted position i i — j, the arches in Ci over C2 and the contribution from a. 
These three terms together are hcj^{i — l,n — k + j — 1) + hc2{i ~ j, k — j ~ 1) + 1. 
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c, ''\\ 




c, '■■■\\ 






i j k 
(a) 


j k i 
(b) 


j ' k 
(c) 



Figure 5. The decomposition used in order to derive the recurrence relation for the 
average height function h{i,n). 



Upon averaging and using ([8]) we obtain the recursion relation for the average height 
function 

^^^^ h{i,n) = ^ h{i,n-k+j-l) + ^ h{i-k+j-l,n-k+j-l) 

i<j<k<n 0<j<k<i 

+ Yl [h{j-l,n-k+j-l) + h{t-j,k-j-l) + l] (14) 

0<j<i<k<n 

In the scahng hmit n —>■ oo, we may insert the scahng ansatz h{i,n) ~ n''T-C{i/n) from 
f fTTj) and replace sums by integrals, which yields after a few manipulations a Volterra-like 
double integral equation for the scaling function 7i{x). It is possible, though tedious, 
to show that the integral equation allows a solution oc — with the scaling 
exponent C = (vTr — 3)/2. Besides the quite comphcated treatment of the integral 
equation we have found evidence for this scaling form from numerical simulations (see 
section [6l) . 

In the sequel, we shall develop a more systematic approach to extract the scaling 
behaviour which is based on recursion relations like (fT4l) . Furthermore, this allows to 
compute sub-leading corrections to the scaling limit and therefore to exactly quantify 
finite-size contributions. 



2.4-2. Generating functions for the local height moments Since the relations for the 
height h are additive in h, it is convenient to deal with the exponential function as a 
generating function of the moments. We thus consider the generating function for the 
height h at site i for a strand of length n 

oo ^ 

E{i, n; z) = (exp {z hc{i, ^))) = XI ^^c(^, ^)^) • (15) 

k=0 

We obtain the recurrence equation for E 

n(n — 1) ^, , \ - 

^ ' E{t,n;z) = E{i,n-k+j - l-z) 

i<j<k 

+ X E{i - k + j ~ l,n - k + j - 1; z) 

j<k<i 
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(16) 



+ e' J2 Eii - l,n - k + J - 1; z)E{i - j,k - J - 1; z) . 

j<i<k 

Note the crucial factorisation in the last term due to the independence of the 
substructures Ci and C2 inside C once the first arch a is chosen. It is convenient to 
introduce the "grand-canonical" generating function 



n=0 i=0 71=0 1=0 

which contains the contribution of strands with arbitrary length n, and which is 
left/right symmetric G{u,v;z) = G{v,u]z). The discrete recursion relation for E 
becomes the non-linear partial differential equation 

+ uv ^ ^ G{u,v;z) 



U 2d' 



V 



2 / g 



u 



u- 



dudv 
1I + 



d 



1 



G{u, V] z) + uv G{u, v; zY 



'It 



[1 — u) \ du J (1— f) \ dv 
We directly derive initial conditions for G{u,v; z) at m = (or f = 0) from the series 
development ( |T7I) . Since by definition the height function vanishes at the ends of the 
strand, we have 



00 00 ^ 

G(0, v;z) = Y, (e^'^cCo.n)^ ^J^^n ^ _ — _ _ 

n=0 n=0 



(19) 



For z = we find G{u, v; 0) = (1 - u)-\l - v)-\ 



2.4-3. Generating functions for the local height h{i, n) and h{n) From G we obtain the 
generating function for the height h{i,n) itself 



d 



OZ 



(20) 



2=0 



F{u,v) = Y,Y.^hc{i,n)), 

71=0 i=0 

Using f|T8l) and G{u,v] 0) = (1 — u)^^{l — f)^\ we conclude that F satisfies the linear 
partial differential equation 



1/2^ 

2 V du"^ dv"^ 



+ uv 



+ 2 



dudv 
u V 
1 — ul — V 



F(u,v) 



u 



d 



u- 



F{u,v) + 



I — u) \ du 

u V 



(1 — f ) \ dv 



9 , 
V— + 1 



(21) 



(l-M)2(l-t;)2 

with initial conditions F{0,v) = F{u,0) = 0. It is straightforward to obtain the 
generating function of the sum of the heights nh{n) (or total area below the height 
curve h{i, n), < i < n) from F{u, v) by setting u = v. 



CO / n 



K{v) := F{v,v) = E Y.(hc{t,n)) v^ = ^^^H^" 

71=0 \i=0 / n=0 

2TI) implies that K is solution of the ordinary differential equation 

2 



1 - vyK"{v) - 2v{l - v)K'{v) - 4(2 - v)K{v) 



(22) 



(23) 
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From h{0) = h{l) = we infer the intial conditions K{0) = K'{0) = 0. Analysis of 
and fl23|) in the hmit u,v ^ 1 will give access to the scaling limits of the height function 
as well as its average h{n). 



3. Mean height and the scaling exponent <^ 

In this section, we derive the exact scaling form of the mean height h{n) from the 
differential equation fl23|) for K{v). In order to get an idea of the scaling limit, let us 
suppose that h{n) scales like 

h{n) ~ cn^ . (24) 



71— +00 



Since ( > 0, insertion of this ansatz into fl22|) implies that the generating function K{v) 
is analytic in the vicinity of f = 0, with convergence radius 1. Its closest singularity is 
situated at f = 1, with a power-like divergence 

Kiv) r. + (25) 

Inserting this ansatz into (1231) . the most singular term is p(C)(l — f)~^~^ = 0, with 
p{C) = + 3C — 2. The roughness exponent ( is thus solution of p{C.) = 0, i.e. 



C. ^ (26) 

We thus identify the roughness exponent with the larger solution 



C = C+ = " ^'^ ^ = 0.561552 .... (27) 

This is the value obtained (using a different argument) by Markus Miiller in [TB]. Let 
us recall that the roughness exponent ( is related to the pairing-probabihty exponent p 
by C + P = 2. Moreover, its inverse is equivalent to the intrinsic fractal dimension (or 
intrinsic Hausdorff dimension) d^ = Thus for our model 



p= L_^l1Z = 1.438447... , d^ = ^ = 1.78077... . (28) 

This exact value for the roughness exponent ( is larger than the one for generic arch 
systems (with weight factors given by the Catalan statistics, i.e. generic trees or branched 
polymers in the dual picture), which corresponds to RNA in the homopolymer phase 
(no disorder), where 

Co = ^, Po = l, 4 = 2. (29) 
However, it is smaller than the value observed in numerical simulations for random RNA 

CandomRNA ~ 0.66 , (30) 

and that of 2-loop RG [12j C = 0.64 for random RNA. 
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We now solve equation (125]) exactly. First, note that a particular solution of the 
full equation is given by 

Mv) = . (31) 

Consequently, we need an appropriate solution Ki{v) of the homogeneous version of 
f l23|) . Performing the transformations 

K,{v) = e-'^{v-lf+'u{z) , z = 2{l-v), C = ^^^^ (32) 

the equation for Ki is changed to a confluent hypergeometric equation for u{z) 

zu"{z) + [2(C + 2)-z\ u\z) - (C + l)u{z) = . (33) 

After a few manipulations, the (appropriate) general solution of this differential equation 
for K{v) = Kq{v) + Ki{v) is of the form 

Kiv) = j^^^, Mi-C, -2 - 2C; 2 - 2v) 

+ C_(l - vY+^MiC + 3, 2C + 4; 2 - 2v) - — ^ (34) 

where M{a,b,z) is the confluent hypergeometric function [18J. The coefficients C+ and 
C_ are fixed by the constraint that K{v) be analytic at f = and that its Taylor 
expansion start at order f^. Hence C+ and C_ are given by complicated and not 
especially enlightening combinations of confiuent hypergeometric functions at z = 2. 
Numerically we find 

C+ = 0.713263... C_ = 0.519299... . (35) 

The first terms of the Taylor expansion of K are rationals 

K{v)=v' + ^v' + lv' + ^v' + ... . (36) 

o o -LO 

The asymptotic limit n ^ oo is equivalent to v —>■ 1~ . In this case, K{v) has a power- 
law divergence and its most singular terms contribute to the leading orders of h{n). A 
Taylor expansion of ( IMI) yields 

c+ , C „ „. C(i-0 „ ..,2 



^-(") - (13^ (^ ^ Ttc - - (1 ; 0(1 + 20 - j - (T^ ^ ■ ■ ■ 

where we have omitted the terms which remains finite v 1^. This expression 
allows to compute the scaling behaviour for the average height h{n) by inversion of 
the transformation. After some algebra, we find for n ^ 1 

Kn) = , , ,, (r(2 + c + n) + c r(i + c + n) 

-1 , 



r(c + 2)nr(n + 1; 

C(i- 0(1 + 



V{C + n) \ - l- n 



(1 + 20 

0.51334 - 1 + 1.31498 n^'^ - n'^ + 0.41413 n^'^ + 0{n^'^) (3^ 
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Therefore at leading order we indeed find the scahng law h{n) ~ cn^ with c = 
C+/r(^ + 2) = 0.51334. . .. Note that in principle, all amplitudes in fl38l) as well as 
subsequent corrections may be computed exactly in terms of hypergeometric functions 
and the gamma function. However, we shall omit these rather lengthy expressions and 
content ourselves with numerical values. This explicit solution will be useful to test 
numerical simulations and the domain of validity for the scaling ansatz (see section [6]) . 



4. Scaling behaviour and scaling functions 
4.I. Scaling form for the F function 

In this section we show that the average-height function h[i,n) = {hc{i,n)) takes the 
following scaling form in the limit of long strands 

{hc{i,n)) = nfrtiix) , x = -. (39) 

This is in fact a particular case of the general scaling form for the moments of h 

{hc{i,nf) = n^'^Hkix) , x=-. (40) 

The partial differential equation fl2Tl) indicates that the generating function F{u,v) has 
in singular lines at m = 1 and at f = 1. These singularities govern the long-strand 
limit n oo with respectively n — i = 0{1) and i = 0{1). The scaling limit n oo, 
i/n = 0{1) is governed by the singularity aX u = v = 1. 

To prove validity of the scaling (l39l) . it is sufficient to show that in the limit M,f — > 1 
the generating function F{u, v) scales as 



-2-C 

with 



F{u,v) = r-2-C^i(cu) (41) 



r = l-^, ^=^, .. = aVr^ (42) 

In terms of the new variables a and r, the scaling limit is r and a —>■ 0, uj = cr^/ = 0{1) 
fixed. In this scaling limit the transformation (120|) (he) F becomes a double Laplace 
transform. The corresponding transformation Tii (x) — » J^i (u) is 

/•oo rl 

J^i(cu)^ dn dxn^+^(l- v^r-r)"^(l + v^r-r)™7^l(x) 
Jo Jo 

/•oo /»1 

^ dn dxn^+^e-(v^+^)"""e-(i-^)™(i-")Hi(x) 
Jo Jo 

pi 

= r(2 + C) / dxHiix) [1- v^+2a;v^]"^'+^^ (43) 
Jo 

To obtain the equation for Ti{uj), we keep the most singular terms in (|2Ti) when u,v ^ 1. 
This gives 



Id Id 



2 \dv? dv^ ) dudv 1 — udu 1 — v dv {1 — u){l — v) 



F{u,v) ^ . (44) 
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Using ansatz fl41l) . we obtain a hypergeometric differential equation for ^1(00) 

uj{l — uj)J^'({uj) + 

Tlius tlie scaling function J-'iiu) is a hypergeometric function 

J-iH = D2Fi(l + C/2,(3 + C)/2,C + 3/2,^). (46) 

where D denotes some constant. We explicitly compute D from the constants C+ and 
C_ obtained in dMD and ([35]) since K{v) = F{v,v) ~ (1 - w)-(2+C)jr^(o) = r-(2+^)j^i(0) 
as f — > 1. One obtains 

D = C+ = 0.713263- ■■ . (47) 

As we shall see in the next sub-section, the form ( l46l) for J-'i{z) implies a very simple 
scaling function for the average height. 



3 + 2C 7 + 2C 



-UJ 



(C + 4) 



(45) 



4-2. Scaling form for the average height function h{i,n) 

Proposition: The scaling limit Tii of the average height distribution h{i,n), as defined 
in (!39|) . is given by a simple "beta-law" with exponent ( 

n^{x) = Ex<{l-x)< , C = ^^^^ (48) 
and the amplitude 

T(2 + 20 

^= a + cm I cr ^^ = 

where C+ is given in fHTj) . 

Discussion and proof: The fact that the average heigth 7ii(a;) scales as x'' for small x 
was already known by Miiller [13]. The simple exact form for 7^i(x) is quite remarkable 
and unexpected. Our first hints for (HH]) came from the numerical simulations that we 
describe in section O 

To prove fHH]) it is simpler to start from (HH]) and to show that it implies the form 
Oin]) for J-'i{z) (the transformation Tii —>■ T\ is linear and one-to-one). Inserting 
into the definition for -F(m, v) when m, f ^ 1 we have (this is equivalent to use (H3l) 

*C 'C poo poo 

F(u,v)^ E V , ^ ^ ^^ u'v^ ~ E / di dj e-*(^-") e-^(^-'') 
+ io Jo 



v/?E r(l-FC)r(2 + C) . F 2 + ^ 2 + 2r-"-^ 

21+2C r(3/2 + C) ^ ^ ^ 2Fi(^l + C,2 + C,2 + 2C,- 



(50) 

i — X/ / 

We now use the quadratic identity for hypergeometric functions [T8| 

2Fi(a,6,26,z)= 1-- (51) 
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in the special case a = 2 + (,b = l + (. We obtain 

v/^Er(i + c)r(2 + c) 



F{u,v) 



21+2? Y{3/2 + C) 



1 - 



U + V 



-2-C 



^ , 2+C 3+C 3 



V — u 



2 — u — V 



14 

(52) 
(53) 



Upon identification with ( l4Ti) (and using the duphcation formula for the F-function, 
[T8] ) we recover the scaling solution (146!) for JF^. Q.E.D. 



4-3. Scaling for higher moments of the height function 

In this section, we study the higher moments of the local height hc{i,n) at site i for a 
strand of length n, {hc{i, nY). Once again, the starting point is a generating function 



Qk 



G,{u,v) = }^}^u^v--'{hc{z,nr) = —G{ 



(54) 



2 = 



EE' 

n=0 i=0 

G{u,v; z) denotes the solution of ffTSl) . From this equation, we are able to recursively 
determine a generating function Gk by a partial differential equation involving functions 
Gk'i k' < k. For example, G2 is a solution of the linear equations 



2^ ^2^ 

2 V 5^2 c}t)2 



g2 ^2 / ^ ^ 

(9M9f (1— m) V ^"^ 



-F(M,t;)+2F(M,t;) 



(1— f) \ dv 



^^7^ + 1 



G2iu,v) 



-G2{u,v) 



(55) 



(l-M)2(l-t;)2 (1_^,)(1_^;) ^ ' ^' V ' ^ ' 

where the right-hand side involves the k = and k = 1 moments. F{u, v) = Gi{u, v) is 
the generating function for the average height {hc{i,n)) studied previously. 



4-3.1. Averaged k-moments Let us first consider the generating function for the 
"integral" of the averaged fc-moment 



00 I. 



K{t-z) = G{t,t-z) = Y,yKk{t) 



fc=0 



(56) 



with 



00 n 



Kuit) = Y.t-Y^{hS,nf) (57) 

ri=0 i=0 

We already know Kait) = (1 — t)^^ and Ki{t) from (IMIl . K{t] z) satisfies the non-linear 
differential equation 

'1 92 t d 2 



2 at2 l-tdt l-t 



K{t,z) =e'K\t;z 



(5^ 



The scaling limit corresponds to t ^ 1 .In this limit, we expect the functions Kk[t) to 
scale as 



Ku{t) - bkil-t)-^-'^' , Ck = kC 



(59) 
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This implies that the average of the fc-moment of the local height scales with the length 
of the strand n as 

- J:ihci^, nf) ^ a, n^' , a, = ^^^^^ (60) 

The coefficients can be computed recursively from the first non-trivial one ai ~ 
0.51334, that we computed previously, see fl38|) . Indeed from fl59l) it follows that K(t, z) 
takes the scaling form 

1 1 °° ?) 

with the scaling variable u = z{l — t)^'^. Corrections are of order (1 — ty~^; considering 
e^K{t, z), they would be of order (1 — t). 

Using ( l58i) and inserting the scaling function IC{u), we obtain up to terms of order 
(1 — t) the equation 

IC{u) + uK,\u) + ^CV/C"(m) = K,{uf . (62) 



We obtain 



+ , ,3 92 + 22^/17 



h2 = K &3 = &? (63) 



Numerically, we find 

/C(m) = 1 + 0.713243m + 0.386756m2 + 0.18727m^ + 0.0851827^^ 

+ 0.0372364m^ + 0.015835'u^ + 0.00659914m^ + 0.00270789m^ + . . . (64) 

As an application, let us evaluate the average height fiuctuations by considering the 
quantity 

^2 = (<^c(fc,^)') - {hc{k,n)f) ^ [a^ - ^^^^±^ ^ 0.055658 n^^ (65) 

We thus conclude that the fiuctuations of the height function remain large in the scaling 
limit (see section [6]). 

4-3.2. General scaling function We now consider the general scaling limit of the 
generating function G{u, v; z) for the moments of hc{i., n). The correct ansatz is 



with 



G{u.,v]z) = T Q{z,uj) (66) 
»i 



Z = ZT ^, T = I , iO = I I (67) 



2 \2-u~v 

In the scaling limit we obtain for Q{z,uj) the equation 

'd^Cod^ , 3-7uj d 1-(1+C)cu d l-Su 
2u;^—-+^S^—+2Cu;t-—+- u;—+ — — —Zt^ + 



duo"^ 2 dz"^ duodz l—uo duo 1—uj dz 1 — uo 

Q{z,ujf (6^ 
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Expanding in z we find the scaling limit for the generating functions of the moments of 
hc{i,n) via a Taylor expansion 



oo ~h 



^(^'^) = Efc!^'^H (69) 

fc=0 



At order k = and = 1 we recover our previous results 

Qoiuj) = , giiuj) = ^i(cu) (70) 

I — UJ 

and for k > 2 recursive second-order linear differential equations for the Qki^) with 
coefficients and second members depending on the previous Qk'{^) {k' < k). From 
Qki^) we obtain the scaling form for the moments 

{hc{i,n)^) ~ n^'^Hkix) , x=- (71) 
The scaling function 7ik{x) is related to Qk{^) by the integral transformation 

gkito) =T{2 + kC) / dxHkix) [1 - v^ + 2xv^"^'+'^^ . (72) 
Jo 

which generalises fj43l) . 

4.4- Simple scaling or multifractality? 

Studying the roughness properties of the height function in the scaling limit, we naturally 
are led to the question of multifractality. We shall argue that within our model the 
height-profile statistics is solely governed by the scaling exponent (. This excludes 
strong fluctuations which might lead to multifractal scaling. To this end, let us consider 
the moments of the local height variations 

Ahk = {\hcit,n)-hcij,n)\') (73) 

In the scaling limit ri ^ oo, we expect a relation of type 

{\hc{t,n)-hc{j,n)\') (X |2-jf'= (74) 

in the regime 1 <^ |i — j] <^ n. U (k = k( there is simple scaling, whereas (k > k( (at 
least for large enough k) implies multifractal behaviour. 

We now argue that we are in the first case and there is no evidence for 
multifractality. Indeed, it is easy to show (using the height picture, and using translation 
invariance to move the point i to the origin of the strand) that the following general 
inequality holds 

{\hc{^,n)-hcU,n)\') < {\hc{i,n)\') , £ = \j - z\ (75) 

We know from (!7T!) that for £ ^ n this scales as 

{\hcii,n)\'') oc n'^Hkii/n) (76) 

In the limit n — > oo, i finite, {\hc{i, n)\^) remains finite, since it is bounded by This 
implies that Hk{x) should behave for small 

Hkix) ~ x''< . (77) 
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This can be shown more rigorously using fl68l) for the generating function Q{z, u) of the 
Qk and the integral relation fl72l) between the Qk and the Tik- The small-x behavior of 
Ti.k{x) is related to the — 1 behavior of Qki^^)- One can check from ( l68i) that the 
function Q{z,uj) must behave when — 1 as 

gCz^co) iM + o(iog(l-c.)) (78) 
Using fl72l) . this implies that 

e^fc(^) ~ ^ Hfc(x) ^ x'^^ (79) 

o)— >i i — a; x^o 

We conclude that 

{\hc{i,n) - hc{3,n)\'') < const. |i - jf^, for K £ < n. (80) 

what implies that ^ ^C- However, we know that Cfc > from general correlation 
inequalities. Hence it follows 

Ck = kC (81) 

what proves the abscence of multifractal behaviour, at least for moments of \hc{i,n) — 
hc{j,n)\. 

4-5. Corrections to scaling 

We can study the corrections to scaling for the height function {hc{i,n)). Let us come 
back to equation (12T1) for the generating function F{u,v) defined by (120|) . A particular 
solution of (12T!) is 



Thus the general solution of ( 1211) is of the form 

F{u,v) = F^{u,v) + C+F+{u,v) + C-F^{u,v) (83) 

where F+{u, v) and F_ {u, v) are two linearly independent solutions of the linear equation 
with no r.h.s. 

q2 ^2 q2 q2 ( d 

+ uv—-- — u— — h 1 



2 du"^ 2 dv'^ dudv {1 — u) \ du 
v"^ f d \ 2uv 



1 — V) \ dv J {l—u){l — v)\ 
It is possible to go to the scaling variable r and u used in equations ( HTl) and ( |66l) 

M = l-r(l + ?/), t; = l-r(l-|/) , ^ = (85) 
and to take for and F_ the solutions which can be written respectively as 

F+{u,v)=t-'-^^:F+{t,u;) , F^{u,v)=r-'~^-T-{T,uj) , C±= ^^~^ (86) 
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(C+ = C is the roughness exponent), such that ^+{t, u) and J-'-ij, uj) have an asymptotic 
expansion in powers of r in the scahng hmit r ^ 0, and are regular in the domain 
UJ G [0, 1[. Indeed (Ell) becomes for J^± the hnear equations 

[-2(r-l)Vc±+2rV(2+C±)-2(r-l)cu(4+2r2+C±+2r(l+C±))] ^±{r,u) 
+2 r-3+2r='(LU-l)'-2C±-2r (cu-1) C±+^ (7+2C±)l oj—^±{t, uo) 
+ [_2r3^2_2 ^ (_2_,_^ _^^) +2(r-l)' (l+r+C±)] ^^-^il^' ^) 

C)2 C)2 

+ (a;-l)4cu2^J-±(r,^)+4(r-l)™^J-±(r,^)+(r-l)V2— J-±(r,a;) =0 (87) 

Note that in the scahng variables the particular solution reads 

F^{u,v) = -T-\l-u)~^ . (88) 

Although not simple, (|87j) implies that its solutions can be expanded in powers of r. 
Indeed, let us expand in r the functions ^±{t,ui) 

°° k ^ 

^±(r,^) = 5^-J-f (^) (89) 

fc=0 



Setting r = in flHTj) fixes the equation for the dominant term to 

2a;(C±+4) J-f (a;)+2 ((7+2C±)a; - (3+2C±)) to ^^±\^) 



+ (c.-l)4c.^^n^(a;)=0 (90) 

This is nothing but the hypergeometric differential equation (H5|) for the scaling function 
jF(u;) obtained previously. Its solution is thus 

T^l\uj) = 2Fi(l + C±/2, (3 + C±)/2, C± + 3/2; ^) (91) 

The expansion in r gives a hierarchy of hypergeometric-like differential equations for 
the correction-to-scaling functions P'^\uj) with a non-zero r.h.s. involving the previous 
scaling functions < k' < k. It is easy to check that these equations admit a 

unique solution T^\uj) which is analytic at = and regular in the domain on G [0, 1[ 
(with a singularity at = 1). 

From this analysis, the coefficients C+ and C_ in the full scaling expansion (!83ll 
are those already calculated in sect. [3], (1351) : 

C+ = 0.713263... and C_ = 0.519299 .. . (92) 

The important result is that, using the inverse transformation F{u,v) {hc{i,n)), the 
average height function takes the general form 

{hc{i,n)) =n^+ (f^n-'^nfix)] + n^- ( n^'' ^L''^(x) ) -1 (93) 

\k=0 J \k=0 / 

where the dominant term is the scaling function obtained in (jlSl) 

nf\x) = Hiix) = Ex^(l - x)^ (94) 
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the leading subdominant term is the last term —1 in fl93p . In fact, since we know that 
the leading order scales like n'^ and therefore only grows relatively slowly with n, the 
correction —1 turns out to be important, even at n = C(IO^), a case that we shall 
consider below. The subleading corrections H^^\x) are distributions on [0, 1] and may 
in principle be computed from the functions T^\uj) via inverse Laplace transforms. 



5. Pairing probabilities 

Up to now we focused on the height function and its scaling laws. However, for the 
original RNA problem pairing probabilities constitute more natural objects. In this 
section, we compute the single-base pairing probability as well as the scaling limit for 
the pairing probability P{i,j). 



5.1. Single-base pairing probability 

As mentioned above, h{l, n) is the probability that a given base is involved in a pair. We 
obtain the generating function g{v) of h{l,n) from F{u,v), introduced in (l20l) . through 
differentiation 

dF{u, v) 



du 



= ^/i(l,n)w" (95) 

«=0 n=0 



According to ( l2Ti) . it is solution of the ordinary differential equation 

(1 - v)g"{v) - 2vg'{v) = g{0) = g'{0) = 0. (96) 

The initial conditions are due to the fact that h{l,0) = h{l, 1) = 0. For our purposes, 
it is sufficient to solve for g'{v) and compare to the derivative of the series expansion 

m- 

g'{v) = -—- = J2nh{l,n)v-' (97) 

^'^ n=l 

Comparison of the series development on both sides leads us to the explicit expression 
Ml, ^) = - E + - E (98) 

k=0 ^ ' fc=0 

In the limit of large strands n — oo, the series converges and we obtain 

lim /i(l, n) = 1 - 4 (99) 

n-^oo e 

For large n ^ 1, the corrections to this result can be determined as follows: 

12 2^ 
h{l,n) = l--- — +r{n) , |r(n)| < — — , (100) 
ne^ n\ mn 

where the bound is obtained by approximating the remaining terms in the sum by an 

integral. We shall reconsider this probability later when comparing the arch deposition 

model to a tree-growth model. 
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5.2. Scaling law for the pairing probability P{i,j) 

In this section we compute the scahng function V as defined in ( ITTI) . First of all, P{i,j) 
is indeed only a function of the distance - despite the fact that we have singled out an 
origin. Therefore we expect for large n the scaling form (fTTj) . 



Pit,j) ~ n-Pv(^-^^] (101) 

n-*oo \ n J 

With this ansatz, relation ^ in terms of the scaling functions for the height field and 
the pairing probability turns into the integral equation 

n^ni{z) = n^-P [ ds [ dtV{t-s) (102) 

This identifies 

p = 2 - C . (103) 

The dimensionless scaling functions for height and pairing probability are then related 
by: 

nz—e pi 

Hiiz) = ds dtV{t-s) (104) 

JQ Jz+e 

Note that we have introduced a small ultraviolet cutoff e in order to circumvent possible 
singularities as t — s ^ 0. Upon differentiating twice, we find the harmless expression 

n';iz) = -iVil-z + e)+Viz + e)) = {V{1 - z) +V{z)) (105) 

It is clear that the limit e ^ will not lead to any problems since potentially divergent 
terms have canceled out. Since V{z) is invariant with respect to the transformation 
z — >• 1 — 2;, it is possible to deduce its exact form 

V{z) = -H'[{z)/2 = E ^ z^-\l - zf-' [1 - C + 2(2C - 1);^(1 - z)] .(106) 

We see that 'P{z) factorises into a beta law with characteristic exponent p = 2 — ( and 
a polynomial correction. Note that a pure beta law is obtained if and only if = 1/2; 
this corresponds to the RNA homopolymer roughness exponent. Since the amplitude E 
is known from poj) . we have entirely characterised the scaling law. 



6. Numerical simulations 



Numerical simulations not only provide a verification of our analytical results in the 
scaling limit n 00 but are useful in order to quantify finite-size corrections. In 
this section, we present a simple algorithm for random generation of hierarchical arch 
structures. Furthermore, we compare the statistics obtained from random sampling to 
the exact solutions. 
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6.1. Outline of the algorithm 

We give a description of the algorithm which we have used to generate hierarchical 
structures C: Information about arches is stored in the "adjacency matrix" ^c{hj)- In 
order to take into account the planarity condition we label each basis i = 1, . . . ,n with 
a "colour" c{i) G Z. During the construction process, two bases i,j may be linked by 
an arch if and only if c{i) = c{j). Furthermore, we introduce two special colours: if a 
structure C contains an arch {i,j) we colour its endpoints with c{i) = 1 and c(j) = — 1 
(what turns out to be convenient). 

The deposition of arches is carried out in the following way: Initially all colours 
are set to c{k) = 0, k = 1, . . . ,n and all entries of the adjacency matrix to ^c{hj) = 

0, z, j = 1, . . . , n. First, we randomly choose a base i among all unpaired bases. Next we 
collect all bases k which may be paired to i without violation of the planarity constraints, 

1. e. with the same colours c{k) = c{i), k ^ i, in an ordered list i. If i is empty, the 
point may be removed from the set of unpaired points and the procedure restarted. 
From the list i of compatible bases we randomly choose a second base j. For simplicity, 
let us suppose that i < j { the converse case is similar). We store information about 
the so-created arch {i,j) by setting ^c{hj) = 1- Moreover we label the starting point 
and the endpoint of the arch with colours c{i) = 1 and c{j) = —1. Finally, in order 
to mark the new substructure due to the insertion of this arch, we set c{k) = i + 1 for 
all i < k < j. We repeat the procedure until no more points can be paired without 
violation of planarity. See figure [6] for illustration of a single cycle. 
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(a) (b) 

Figure 6. A single step of tlie construction of liicrarchical structures: (a) random 
choice of the point i. Point j is chosen amongst all points of the same colour as i (light 
gray), (b) Once the arch is determined all the points i < k < j are re-coloured 
in order to mark the new substructure (dark gray). 

Once this procedure is finished, the matrix ^c{hj) contains all information about 
the structure. To compute the height field, the colours c{i) may be used as well: if 
we set c{i) = for all i such that c(i) > 1, then we obtain a sequence {c(i)}"^^ with 
entries 0,±1. It precisely corresponds to the discrete derivative of the height function 
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c{i) = hc{i,n) — hc{i — ^,n). Therefore, we can reconstruct 

i 

hc{i,n) = Y,c{k) (107) 
fc=i 

For a given strand of length n, we perform this construction times in order to average 
over the samples. This algorithm is a variant of the point process to be discussed 
in section [71 Though not being dynamically equivalent to the arch deposition model, 
the key feature of partitioning into independent sub-structures leads to the same final 
probability law in configuration space. 

6.2. Results 

We have constructed structures with up to n = 6500 bases in order to test our theoretical 
predictions. For n < 200 bases, we have sampled 10^ structures whereas for n > 200 
bases, 10^ structures per data point were sampled. 

Single-base pairing probability h{l, n) : Results for the probability that a base is paired, 
which equals the height h{l,n), is presented in figure [7^. We find agreement with the 
theoretical prediction from (|98|1 within errorbars. 




(a) (b) 

Figure 7. (a) Deviation of the single base pairing probability h{l,n) from 1 — e^^ as 
a function of n. The dashed line is the theoretical prediction (|98p (b) Log-log-plot of 
the average mean height h{n) as a function of the number of bases for 5 < n < 6500. 
The straight line corresponds to the theoretical prediction from (|38p . the dashed line 
indicates the scaling limit. (The error bar is of point size.) 



Averaged mean height h{n), and its fluctuations: In figure [7b we compare results for 
the averaged mean height h{n) to the theoretical prediction in the limit n oo. Taking 
into account all terms of ( l38l) is sufficient to show that the difference between numerical 
results and theory is of the order of the statistical error, see figure [8^. 
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(a) (b) 

Figure 8. (a) Deviations of the averaged mean liciglit h{n) from tlie tlieoretical 
prediction, normalised by the result of . Clearly, the deviations are of smaller than 
the errorbars. (b) Log- log-plot of the averaged second moment A2. (The error bar is 
of point size.) The straight line presents the scaling limit. 



In order to compare our data to the results of section l4.3.1l on averaged fc-moments, 
we evaluate the height fluctuations via A2 = Ylk=ii{^c{k, nY) — {hc{k, n))'^)/n in flgure 
[8]d. The data shows good agreement with the prediction from (!65l) . 

Averaged height function h{k,n). Results for the averaged height functions are shown 
in flgure [9H,b. In order to point out universal behaviour we plot h{k, n)/n'' as a function 
of X = k/n. In flgure [9)d, we compare the data to the flrst-order corrected scaling limit 
n~^l-L{k/n) — 1 where 7i(x) denotes the scaling function from fHHj) . The deviations 
Ah{k, n) = h{k, n) — n~^7i{k j n) + 1 are large at the ends k = 1 and k = n. 




Figure 9. (a) Scaling plot for the height function 7Y(fc, n) for two strand lengths 
n = 100 and n = 6000. The dashed lines correspond to the scaling function plus the 
first finite-size corrections to the scaling limit n"^ (n^7i(fc/n) — 1). (b) Residual scaled 
deviations Ah{k, n)/n'^ for different strand lengths n — 50, 100, 200. 



We examine the deviation of the height function h{k, n) from the scaling limit 
by evaluation of its value at k = n/2. Figure [10] shows the rescaled deviation 
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n'''^{h{n/2,n) — n'^l-i{k/n) + 1). Numerically, we find 

n-%{n/2,n) - n^n{k/n) + 1] = 0{n-^] 
in agreement with the scaling form (!93|) . 
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Figure 10. Log-log-plot of deviations of the height function h{k,n) at fc = n/2 as a 
function of n. The straight line corresponds to the function to n~^. 



7. Growth model 

The model considered in the previous sections (model A) is a deposition model. The 
size of the system is fixed; to study the folding of a strand with n bases, we start from 
a set of n unoccupied points {!,■■■ ,n} on the line and successively deposit arches in a 
planar way until the system is full (no deposition possible). Systems with different size 
n and n' are a priori different. 

In sections 17.11 and 17. 2^ we show that this arch deposition model A is equivalent 
to a stochastic growth model G for arch systems, where we start from a system with 
no points. At each time step t we deposit a new point according to a simple stochastic 
process, and create a new arch whenever it is possible. We show that in the growth 
model G the statistics for the arches at time t is the same as the statistics of arches of 
the deposition model A for a system of N = t points. 

In section 17.31 we shall also show that this stochastic growth process can be 
reformulated (by a simple geometric duality) as a tree growth process T. In section [731 
as an application, we compute in a simple way local observables of these models, such as 
the asymptotic (at large time) distribution of the number of branches for a vertex of the 
growing tree, which is related to the asymptotic distribution of substructures (maximal 
arches) in the arch model. Finally, in section I7.5[ we study the dynamics of this growth 
model, and compute time-dependent pairing correlation functions. 

7.1. Arch growth models via point deposition 

7.1.1. Closed-strand growth model G Let us first define the model G (closed model), 
illustrated on figure [ID 
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At time t = we start from a closed line (a circle) with no point. At time t = 1 we 
deposit a point on the circle. Assume that at time t we have already deposited t points 
and constructed a maximal planar arch system between these points. Namely there are 
riait) arches and nf(t) = t — 2na(t) free points such that it is impossible to construct a 
new arch linking 2 free points without crossing an already constructed arch (planarity 
condition) . 

At time t + 1 we deposit a (t+ 1)*'^ point, with equiprobability 1/t on the t intervals 
separating the t already deposited points. If it is possible to draw a planar arch between 
this last point and one of the free points (i.e. an arch which does not intersect one of 
the existing arches) we add this arch (It is clear that this arch is unique, otherwise the 
existing planar arch system at time t would not be maximal). Otherwise the new point 
stays free. 




Figure 11. Model G: successive deposition of points on a circle. Unlinked points are 
marked in white, linked points in black. 



Here we give all possible configurations from n = 2 up to n = 9 vertices, with their 
pr obabilities. Furthermore figure [12] shows a sample for n = 600 points. 
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(114) 



Figure 12. An explicit example of a structure with n = 600 points. 



7.1.2. Open-strand growth model G' Let us now define a slightly different model G' 
(open model). 

At time t = we start from an open line with no point. At time t = 1 we deposit 
a point on the line. At time t, we assume that we have already deposited t points on 
the line and constructed a maximal planar arch system between these points. Namely 
there are na{t) arches and n/(t) = t — 2na{t) free points such that any link between two 
of these free points necessarily intersect one of the na{t) existing arches. At time t + 1 
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we then deposit a. (t + 1)*^ point, with equiprobabihty l/(t + 1) on the t + 1 intervals 
separated by the t aheady deposited points. If it is possible to draw a planar arch 
between this last point and one of the free points we add this arch, otherwise the new 
point stays free. 

Note that model G' can also be viewed as model G with an additional inactive 
point, marking the cut. 




Figure 13. Model G': successive deposition of points on an open line. Unlinked points 
are marked in white, linked points in black. 



7.1.3. Relation between model G and model G' It is clear that a configuration C of 
model G can be obtained from a configuration C of model G' by closing the line and 
that all configurations C that are equivalent by a discrete rotation give the same C (see 
figure [T4l) . In other words, the configurations C of model G are the Z„-orbits of the 
configuration space of G' under the action of discrete rotations. 




Figure 14. From a configuration C of G' to a configuration C of G 



7.2. Equivalence between the growth model G' and the deposition model A 

It is clear that the arch configurations C of model G' are the same as the arch 
configurations of model A. It is less obvious that the probability for each configuration 
in both models are the same. 

Theorem: The probability P{C) of any configuration (i.e. class of diagrams) C in models 
G' and A are the same. 

Pa{C) = Pg'{C) . (115) 

To prove the theorem we start from the recursion equation ([8]) for the configuration 
probabilities in the arch-deposition model A, that we obtained in section [2m and rewrite 
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here for completeness 

^^(^)= E -r^, Pa{Ci) Pa{C2) (116) 
n(n — 1 

arch aeC ^ ' 

This recursion relation, together with the initial condition P = 1 for the n = and 
the n = 1 configurations (no point and a single free point), is sufficient to obtain all 
probabilities. 




Figure 15. Probability recursion (|116p as decomposition of a configuration C in model 
A 



We now prove that the probabilities in model G' obey the same recursion relation. 
For this we first need to relate the propabilities in model G to those in model G'. 

Lemma: Let C be a configuration with n points in model G' (successive deposition 
of points on a line) and C its equivalent configuration in model G (successive point 
depositions on a circle). Let s{C) be the symmetry factor of the configuration C, i.e. the 
number of cyclic rotations that leave C invariant. Then 

^G'(C) = ^^g(C) (117) 

Proof of the lemma: It is clear that in model G' any deposition process of n points 
on the line is uniquely specified by the bijection i x{i) where x{i) is the position at 
time t = n of the point deposited at time i. x is a bijection on {1, n}, i.e. a permutation. 
Any process is equiprobable, therefore the probability for any x is p{x) = l/n\. It is 
equivalent to successively create the arches as soon as this is possible, or to create all 
the arches at time n, with the constraint that any point x{i) can only be connected to 
the points x(j) with j < i. To any permutation x is associated a unique arch system C 
and the probability for C is 

Pg'(C) = — number of x ^ C = — cardix : x C} . (US) 
nl nl 

Two configurations C and V of model G' are equivalent in model G if they are equivalent 

by some Z„ rotation r. 

C = V C = V. (119) 
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This means that if x is a permutation for C, ?/ = r o a; is a permutation for V. Hence 
there are as many permutations for C as for V. 

C = V Pg'{C) = PaiV) . (120) 

We now count the number of C which are equivalent by rotation and give C. This is 
obviously 

11 

number of C — >■ C = — — . (121) 

Now we go back to the model G. Any point deposition process on the circle is also 
in bijection with a permutation, but now with one point fixed, for instance x(l) = 1. 
Therefore 

Po{C) = J2Pa'{C) = -^Pg'{C) . (122) 




Figure 16. The decomposition of a configuration of model G used in the proof. To 
be compared with figure [T5l 



We now go back to the proof of the theorem. In model G any configuration C (with 
n points) can be constructed by first depositing a couple of points (1,2), which form 
a first arch ai and then by depositing rii points to the right of a and n2 points to the 
left, with of course ni + n2 = n — 2. Let us denote Ci and C2 the arch configurations 
to the right and to the left of oi in C. These configurations are arch configurations of 
model G', not of model G, since the arch oi cuts the circle into two segments. Once the 
first two points are deposited, amongst the (n — 1)! possible ways to deposit successively 
the last n — 2 points, each either to the left or to the right of Oi, there are (n — 2)! 
possible ways to deposit ni points to the right and n2 points to the left, independently 
of (ni,n2). In other words, the distribution for (ni,n2) is uniform. 

prob(ni,r;,2) = — ^— , Ui + n2 = n - 2 . (123) 
n — 1 

This can be shown easily by using the recursion relation 

prob(ni, 712) = prob(ni - 1, 77,2) — + prob(ni, 77,2 - 1) — , (124) 

'^i + ^2 + 1 ni + n2 + l 
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with initial condition prob(0,0) = 1. Once this is done, the conditional probabilities to 
obtain Ci and C2 are independent, and given by Pg'(Ci) and Pg'{C2)- 

The total probability to obtain a configuration C in model G is therefore given by 
a sum over all (first) arches ai in C. Each term of the sum is the probability that ai 
is the first deposited arch, and that one obtains Ci and C2 in process G'. There is a 
counting factor 2/s(C) associated to each initial arch oi, where the factor of 2 accounts 
for the two possible choices for the first point 1 on oi, and the symmetry factor l/s{C) 
is there to avoid multiple counting when several arches are equivalent. Therefore we 
have finally 



Using Lemma flll7p the symmetry factor disappears and we obtain for the probability 
in model G' the recurrence equation 



This is exactly the same recurrence relation as for model A. The initial condition are 
the same for n = and n = 1, which proves the theorem. 

7.3. Equivalent tree growth processes 

7.3.1. Duality with trees There is a well-known dual description of planar arch systems 
in terms of planar trees. Represent faces by vertices, and arches by links between two 
vertices. In our model, we have also free points deposited on the external circle but not 
yet linked to another point by an arch. Every face of the planar arch system has at most 
one such free point. We represent such a face by a white vertex o with a white arrow 
pointing towards the free point. Every face with no free vertex is represented by a black 
vertex •. We thus obtain a dual description in term of decorated planar trees with at 
most one arrow per vertex (see figure [T7|) . Within this dual description, the model G is 
a planar tree growth model T defined as follows. 

7.3.2. Tree growth processes T At t = we start from the tree with a single black 
vertex and no link. We define the tree-growth process as follows: As illustrated in figure 
[T8l at each time step we 

- either add an arrow to any black vertex, so that it becomes a white vertex (for a black 

vertex with k links, i.e. a fc- vertex, there are k different ways to add an arrow); 

- or add a second arrow to a white vertex (for a white vertex with k legs there are k + 1 

different possibilities); and then transform this white vertex onto two black vertex 
with a new link ortogonal to the two arrows. 

See figure dS] for an illustration. This internal budding process is a specific feature of 
our growth model. 




arches aSC 



(125) 




(126) 
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Figure 18. Elementary growth steps for the decorated tree model: at each step we 
add an arrow to some vertex; if there is already an arrow the vertex splits in two. 
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7.3.3. Another tree growth process T' A similar growth process is obtained if we forget 
about the arrow position for o-vertices. One considers trees with black vertices •, and 
white vertices o if there is an arrow (see figure [T9l) . Indeed it is easy to see that at 
each step the position of the arrow around a white vertex is equiprobable, i.e. there is a 
probability 1/k for an arrow to be at a given position on a type o fc- vertex. With this 




Figure 19. The position of a arrow around a o vertex is uniformly distributed. 

property, we consider undecorated trees made out of •- and o-vertices. We start from a 
single • vertex at time t = 0. At each time step we 

- either transform a black A;-vertex into a white vertex with probability weight w,^o = k 

(where k is the coordination number of the black vertex); 

- or transform a white /c-vertex into a pair of black vertices, one fci-vertex and one k2- 

vertex, with ki + k2 = k + 2 (and ki and > 0), with a uniform weight Wo^,» = 2/k 
for each occurrence (since for each ordered pair {ki, k2) there are k/2 possible ways 
to split o — >• • • and there are k + 1 such ordered pairs) . 




Figure 20. Growth processes for a /c = 4 vertex in the undecorated vertex model. For 
k = 4 the process • ^ o has probability weight w{k) — k — i, each process o ^ • • 
has w{k) = 2/k = 1/2. 

2 

w,^o(k) = k and Wo~,,,(k) = — . (127) 

k 

Let us denote the number of •- and o-vertices at time t by n,(t) and Uoit) respectively. 
For any tree created through this growth process up to time t, it holds the Euler relation 

2n.{t) + 3no{t) = t + 2 , (128) 
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This relation proven by induction. For t = 0, n, = and Uo = 1. During a time step 
t t + 1 we either have (n,, no) —>■ {n, — 1, rio + 1), or (ra,, Uo) {n, + 2, rio — !)• In 
both cases 2n, + Suo increases by 1, as does t + 1. 

The transition probabihty pc^c to go from C C" is defined from the probabihty 
weights wc->c' by 

Pc^c' = ^ • (129) 

2^(^„ wc-^c" 

Thus the probabihty P{C, t) at time t to be in a configuration C is obtained recursively 
by 

P{C, t) = Tpc-^c P{C\ t - 1) = V J^^'^^ P{C\ t-l). (130) 

Note also that the dual of an arch configuration C in model G' is a rooted tree 
footed (decorated with white arrows, or with black and white vertices as explained 
above). Thus model G' is dual to a growth model for rooted trees (see figure [2T|) . 





Figure 21. Open planar arch systems are dual to rooted trees. The root is indicated 
by the barred line pointing to the top. 



7.^. Mean-field calculation of local observables 

In this section we present a mean-field theory for the point deposition model G. It 
allows to compute local observables defined on the tree structures such as average vertex 
densities, coordination numbers etc. in the large-size or long-time limit. More precisely, 
the mean-field approximation amounts to neglect fiuctuations in the limit t — ^ oo which 
are of order 1/t. This allows to transform the problem into a stationary process. 

7.4- i- Vertex densities The simplest observables are local observables, such as the 
average number of vertices of a given type. Let us denote by n,{k, T) and rioik, T) the 
total number of fc-vertices of type • and o in a tree configuration T obtained at time t, 
starting from • at time t = 0. We have shown above in (I128P that at any time t and for 
any tree configuration T 

oo 

^(2n.(fc,r) + 3no(fc,r)) =t + 2 . (131) 

k=l 

Starting from T at time t, the total number of weighted moves t — >■ t + 1 (i.e. of ways 
to add a new point on the dual configuration) is 

oo 

S(t) = {kn,{k, r) + {k + 1) no{k, T)) = t . (132) 

k=l 
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For simplicity we denote by bk{t) = {n,{k,T))^ and Wkit) = {no{k,T))^ the average 
number of black and white /c-vertices at time t. We write down a master equation for 
their time evolution during a step t ^ t + 1. To this end, we have to evaluate the 
transition probabilities for transformations of black and white vertices. 

During the time step t ^ t + 1 the probability for a given black /c- vertex to become 
white is 

Similarly, the probability for a given white fc-vertex to split into a pair of black ki- and 
fc2-vertices (with ki + k2 = k + 2) is 

p{°k •fci'te) = = ■ (134) 

Hence the master equations for the vertex numbers are given by 



hit + 1) =bkit) + -\-khit) + 2 Wgit)] (135) 

q>k 



U-kb,{t)+2 J2 

\ q>k~l 

Wk{t + 1) = Wk{t) + -^{k hk{t) -{k + l) Wkit)) . (136) 

In the large time limit t — > oo we expect the vertex numbers bk{t) and Wk{t) to be 
extensive, i.e. proportional to t. Therefore, we define (assuming that the limit exists) 
the density of black and white vertices as 

(3k = liiii ^^-^ and ujk = hm ^ . (137) 

t^oo t t—tOO t 

Consequently, from the master equations we find two coupled recurrence equations 

2 . ^ k 

(3k = ^— Y 2^ , ^k = -^-^ (3k (138) 

q<k-l 



whereas relation f ll32p implies 

Y^kPk + {k + l)i^k = l . (139) 



The solution of these equations for the densities is (see Appendix Appendix A for the 
derivation) 

12^ 1 k 2^ 

= 7-Jk^v = ^ (^ • 

7.Ji..2. Results for vertices and related local observables The exphcit expressions for the 

vertex densities allow to determine some interesting quantities. For a system with t 

bases, the average number of black and white vertices are 

g2 _ 3 

n,{t) = tx22f^k = tx and (141) 

k>0 

no{t) = t x^ujk = t X ^ . (142) 

k>0 ^ 
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Therefore, the average number of vertices is given by 

n,{t) + no{t) = t X [1 - e-^) (143) 

We are aheady famihar with the expression on the left-hand side: because of the duahty 
between trees and arch diagrams, we have just calculated twice the average number of 
arches in the large strand limit. However, this is nothing but the number of bases (here 
t) times the single base probability (1991) . In fact, this observation is consistent with the 
value for the fraction of white vertices 

uj = — = ^ujk = e-^ = 0.135335 (144) 

fc>0 

because of the relation lim„_,oo h{l, n) = 1 — u. 

In order to learn more about the average tree structure, we compute the average 
coordination numbers in mean-field theory: 

E^A 2, 1 E^^fc 2_o 

{k.) = '-^^ = ^ = 1.91136, ih) = '-^ = ^ = 2.19453 . (145) 

k>0 k>0 

On average vertices have two legs. The probability for a branching, i.e. the probability 
to have a vertex with at least 3 points is 

^h^'t'l' = = 0.269584 . ( 146) 

}2k=i^k + f3k 3(e2-l) 

More specifically, the probabilities p{k) to have a branching with a black or white k- 
vertex are: p(l) = 0.41738, p{2) = 0.313035, p(3) = 0.166952, p(4) = 0.0695634, 
p(5) = 0.0238503, .... We thus conclude that branchings (i.e. vertices with at least 
three legs) are not rare. 



7.4-3. Substructures and exterior arch statistics The arch-tree duality allows us to use 
the tree growth model to analyse the number of substructures of arch diagrams. A 
substructure is defined as a maximal (or exterior) arch which has no further arch above 
itself (see figure 1221) . 




Figure 22. Example of a configuration with k = 3 substructures and white-coloured 
root. 
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We characterise an arch diagram by {k, a) where k denotes the number of substructures 
(number of maximal arches) and cr = • or o if the root vertex is black or white. We are 
interested in the large-time probability distribution p{k, a) that the "state" of the arch 
system is {k, cr). Consequently, the probability that the arch diagram has k substructures 
is given by p{k) = p{k, •) + p{k, o). A {k, •) state is dual to a tree with a black k- 
vertex with a marked leg. The same holds for {k, o) states. Taking into account the 
combinatorial factor of k for marking a black fc-vertex (or equivalently cutting open a 
circle at a /c- vertex), and a similar factor of + 1 for a white circle (remembering that 
the additional white point, see figure [221 allows for one more option to cut the circle 
open) the probabilities p{k, •) and p{k, o) are proportional to the fraction of black or 
white vertices respectively 

l^ki'^Pk + {k + l)ujk) 

P{k, o) = ^..^^tr^-u'n ^ = (^ + 1)^^ ' (148) 

l^kV^Pk + [k + l)ujk) 

where we have used the Euler relation ^^(/K/Sfc + {k + l)ujk) = 1 to simplify the results. 
Using (11401) we obtain 

2^ 1^ 
= e^(fc + l)(fc-l)! °) = e^(fc + 2)(fc-l)! - ^'''^ 

With the probability distribution for the number k of substructures 

p{k) = pik, .) + pik, o) = , (150) 

we are able to evaluate its moments (k"^) = k"^p{k) in order to characterise the arch 
diagrams. The average number of substructures is {k) = (5e^ + l)/(2e^) ~ 2.56767 to 
be compared with the result (^) Catalan ~ found for Catalan structures [19]. For its 
variance we find (A;^) — (A;)^ = (9e^ — 16e^ — l)/(4e^) ~ 1.70408 which is smaller than 
the corresponding value (^^) Catalan ~ (^'^ Catalan ~ ^' therefore conclude that 
hierarchically constructed structures fluctuate less than generic, equiprobable Catalan 
structures. 



7.5. Dynamical correlations 

We can also compute dynamical quantities in the tree growth model G (and G'). The 
dynamics of this model is interesting in its own, but note that its dynamics is different 
from the arch deposition dynamics of model A. Let us give a few examples. 



7.5.1. Probability of non-immediate pairing Consider process G. Having deposited a 
point at time t, one might ask for the probability that it does not get paired immediately 
with a free point already present. This equals the probability that we add an arrow to 
a black vertex, not to a white one. At large t it is 

^ A; = = 0.567668 . (151) 
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7.5.2. Time- dependent pairing propabilities What is in model G the probabihty "^{i,]) 
that the point deposited at time ti = i is paired with the point deposited at time 
t2 = j > i, as a function of i and j? This amounts to the following event. At time 
t = ti = i a point is deposited on the circle so that no arch is formed, i.e. a certain black 
fc-node is converted to a white k-node. The probability of this event reads bk{i) k/i, see 
(11331) . This particular node then remains white up to time t = t2 = j where it is 
converted to a pair of black nodes. In the timestep t —>■ t + 1 < j the probability of 
keeping the white k-node unchanged is 1 — {k + l)/t, the probability of splitting it is 
{k + l)/t, see (11341) . Thus, if we start from a fc-node at time t = 1, the probability is 



hk{i) k 



n 

.t=i+l 



1 



k + l 
t 



k + 1 
3 



and the total probability is 



hk{i) k 



n 

t=i+l 



1 - 



k+l 
t 



k + l 
3 



It is interesting to consider the large-time, i.e. large-size limit t 
i/j = 0{1). Indeed, using bkit) ~ t/^fc and 

i-l / , 1 \ / - X k+l 



(152) 
(153) 

oo with i,j —y oo. 



n 

t=i+l 



k+ 1 

t 



the time-dependent pairing probability takes a simple scaling form 



oo ^ 

.k=l 



k+l 



{k 



g2 j3 



J 



(154) 



(155) 



Note also that in the large-time, i.e. large-size limit t 
time i gets paired with probabihty one. Indeed 

l + e2 



2e2 



0.567668 



oo a point deposited at a finite 



(156) 



is the probability (I15ip of non-immediate pairing. 



8. Conclusions and outlook 



To summarise, inspired from the subject of RNA folding, we have introduced and studied 
a growth model of planar arch structures (which can be viewed as an arch deposition 
process) . The construction of arch structures is similar to processes generated by greedy 
algorithms. The arch-growth model turns out to be amenable to analytical calculations. 
We have calculated the generating functions for the local height, and their moments. 
This allowed us to obtain the scaling exponent ( for the height, the exponent p for 
the pairing probability, the corresponding scaling functions in the limit of long strands 
72 — i> oo, as well as finite-size corrections. We also proved the absence of multicriticality. 
These results were then confirmed by numerical simulations for systems of sizes up to 
n = 6500. 
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In a second step, we have defined an equivalent tree-growth model. This model 
involves growth by vertex splitting as well as by vertex attachment. This growth process 
allows to generate RNA configurations with arbitrarily large strands (number of bases). 
This allows us to obtain quantities as e.g. the probability, that a point gets paired, 
analytically. 

This work leaves open many interesting questions: 

- Some properties (e.g. distances on the tree, fractal dimension) are easy to study 
in the arch-deposition formulation, while some other properties (e.g. substructure 
statistics) are easier in the tree-growth formulation. It would be interesting to have 
a better understanding of this fact. 

- The equivalence between the arch-deposition process and the tree-growth process 
is very specific to models A and G. We have not been able to find a tree-growth 
process which is equivalent to the compact arch-deposition model A, although this 
model A is in the same universality class as the non-compact arch-deposition model 
A. 

- Is there a tree-growth process which gives the statistics of planar arches in the high 
temperature phase where all arch structures have the same probability (i.e. the 
statistics of the so-called "generic trees" or mean-field branched polymers)? 

- Arch structures and trees appear in many problems in physics, mathematical 
physics, combinatorics, computer sciences, etc., in particular in integrable systems 
(Razumov-Stroganof conjecture, loops models), random permutations, random 
matrix models, and interface growth. Are the kind of models introduced in this 
article related to these problems? 

Finally, since our scaling exponent ( = (vT7 — 3)/2 deviates from the value found 
for random RNA ( 0.66, we conclude that the low-temperature phase of random RNA 
is governed by rules which are more complicated than the greedy algorithm. It would 
be interesting to find a refined scheme that yields statistics closer to random RNA in 
order to comprehend the nature of the glassy phase of random RNA. 

Acknowledgements 

This work is supported by the EU ENRAGE network (MRTN-CT-2004-005616) and the 
Agence Nationale de la Recherche (ANR-05-BLAN-0029-01 & ANR-05-BLAN-0099-01). 
The authors thank the KITP (NSF PHY99-07949) , where this work was started, for its 
hospitality. We are very grateful to M. Miiller for stimulating discussions, and providing 
us with a copy of his PhD thesis. We also thank R. Bundschuh, P. Di Francesco, T. 
Jonsson, L. Tang and A. Rosso for useful discussions. 



A growth model for RNA secondary structures 39 
Appendix A. Mean field equation for the vertex densities 

This appendix contains a detailed presentation of the computation of the vertex densities 
within the framework of mean field theory. We start from the Euler relation 

oo 

J2{kf3k + {k + l)ujk) = l (A.l) 

k=l 

and 

q=k — l 

Taking these two last relations together, we find 

2 = '''-'^ 2 = k+1 + 2 ^^-^^ 

^{k + 2fPk+i = {k + 2){k + 3)Pk+i + 2kpk (A.4) 

This leads to the following differential equation for the generating function B{z) = 

(^^^ + 1^ {B{z) - (3^z) = A (^^^ + 1^ (^(^) _ p^, _ f3^,2^ + 2z'B'{z) (A.5) 

Using the fact that P2 = 2/3i/3, one finds 

z{z - l)B"{z) - {2z^ -3z + 2)B'{z) + B{z) = -2/3i (A.6) 
whose solution is 

Biz) = -2/3, + C,— + C2 . (A.7) 

The simple pole at 2; = is removed via setting C = Ci = — C2. Furthermore, 
consistency requires B{z = 0) = what yields C = Pi/2. Thus, the generating function 
is determined up to a factor: 

The overall factor (3i is obtained by insertion of these expression into lA.li The result is 
so that /3i = e~^. Therefore, we obtain the vertex densities 
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Appendix B. The compact arch deposition model A 



In the compact arch deposition model one deals with strands with an even number of 
bases i = 2n, and the arches are always between an even and an odd base a = (even, odd) 
or (odd, even). At the end of the deposition process, there are no free base and there 
are always n = i/2 arches. The recursion relation ([8]) for the probabilities -Pa(^) the 
configurations C becomes in this model 



Z — / 

arch aSC 

The recursion relation for the generating function of the height 



(B.l) 



oo e=2n 

E E 

n=0 i=0 



(B.2) 



is easily derived and reads 
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(B.3) 



(to be compared with the equation ( l2Ti) obtained for the non-compact model A). 

The scaling limit £ — oo is still given by the singularity at M,f — >^ 1. In this limit 
the dominant (most singular) terms are 
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F^{u,v) (B.4) 



2{1 -u)du 2{l-v)dv {l-u){l-v) 
This is the same equation than equation (144|) for the scaling limit of the generating 
function F{u,v) for the non-compact model A. Therefore the scaling limit for the non- 
compact model A and the compact model A are the same. The same result holds for 
the higher moments correlation functions and the A^-points correlators. 

Finally, let us mention that, although the deposition models A and A are very 
similar, we have not been able to construct a growth model (i.e. a point deposition 
model) which could be equivalent to the compact arch deposition model A. 



Appendix C. Multicorrelators 

We can extend the recurrence equations (IT^ and (ITSl) to compute correlation functions 
for heights at several points of the strand. Let us consider the 2-point correlators. They 
are the expectation values, at two points i and j, of 

{hSi,n)''^hcit2,n)''') . (C.l) 

A generating function for these correlators is 

G2{u,v,w;zi,Z2) = Yl u't;^-'u;"-^'(e"^'^^(''")e"^^'^(^'")) (C.2) 

0<j<j<71 
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Figure CI. The six different ways to deposit the first arch (fe, I) w.r.t. the two points 
i and j correspond to the six different terms in the r.h.s. of (|C.4p 



The recurrence equation is obtained by considering all the possible positions for the 
first arch (/c, /) with respect to the two points i and j (see figure [CT]) Denoting by C the 
length operator 



„ I d d d 

ou ov ow 



(C.3) 



and G{u, v; z) the 1-point function studied in Sect. 14. 3[ we obtain the linear PDE for 

G2 



-C{C-l)G2{u, V, w; zi, Z2)-- 



u 



1 — n V du 



d \ V 



1 — f \ dv 



d \ w 



9 , 
w- — hi 



+Mf e^^ G{u,v; zi)+vwe^'^ G{v,w; Z2)+uw e^^^^^ G{u,w; zi + Z2) 



1 — w \ dw 
G2{u,v,w;zi,Z2){CA) 



Each term in the r.h.s. of (10.4^ corresponds to one of the positions in figure ICll (in 
the same order). The boundary conditions are given by the cases m = 0, f = or 
w = where the 2-point function G2 reduces to a 1-point function G = Gi or a 0-point 
function Go- 

Similarly, we may consider the 3-point function 

G3{u,v,w,x;zi,Z2,Z3) = J2 (mV-%'^-'x"-V^^(*'")+"2^(^'")+"^'^('='")) 

0<i<j<k<n 

G3 satisfies a linear PDE with coefficients involving the 1-point and 2-point functions 
G and G2. 
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